function rad2brt_f(np, wavenumbers, radiance) result(brt)
    use iasing_osse
    implicit none
    ! Arguments
    integer(JPI4), intent(in)   :: np
    real(JPR4), dimension(np), intent(in)   :: wavenumbers, radiance
    real(JPR4), dimension(np)   :: brt

    ! Auxiliar variables
    real(JPR8), dimension(np)   :: rad8, brt8
    real(JPR8), dimension(np)   :: wnu8_scaled
    real(JPR8), parameter       :: scale = 100

    wnu8_scaled = real(wavenumbers, 8)*scale
    rad8        = real(radiance, 8)

    if (rad8(1).gt.1) then
        rad8 = rad8/1000
        !print*, "Process: rad2brt_f"
        !print*, "Correcting scale factor because bad radiance units..."
        !print*
    endif

    brt8 = CST_HCK * wnu8_scaled / &
           log((scale * CST_2HC2 * wnu8_scaled**3 / rad8) + 1)
    brt = real(brt8, 4)

end function rad2brt_f
